function Verify_Solution_Model_20220926(parms,gesol)

    NN = parms.NN; Nx = parms.Nx; Nz = parms.Nz;
    CV = zeros(NN,Nx,Nz);
    S = NaN*zeros(NN,Nx,Nz);
    
    Nlb = parms.grid_N(1); Nub = parms.grid_N(end);
    xlb = parms.grid_x(1); xub = parms.grid_x(end);
    
    [N_mesh,x_mesh] = meshgrid(parms.grid_N(:),parms.grid_x(:));
    
    zeta_fun = @(N,phi)(N + (1-N).*(phi.^(1-1/parms.chi))).*((N + phi.*(1-N)).^(parms.gamma-1));
    phi_of_x = @(x)1./(1+exp(-x));
    
    C_next = zeros(Nz,1);
    S_next = zeros(Nz,1);
    x_next = zeros(Nz,1);
    phi_next = zeros(Nz,1);
    SDF_full = zeros(NN,Nx,Nz,Nz);
    SDF = zeros(Nz,1);
    
    for iz = 1:Nz
        expected_znext = parms.StateTransitionProbs(iz,:)*parms.grid_z(:);
        if parms.OPTION_x_standardize_shock
            std_znext = sqrt(parms.StateTransitionProbs(iz,:)*((parms.grid_z(:)-expected_znext).^2));
        else
            std_znext = 1;
        end
        for ix = 1:Nx
            for iN = 1:NN
                
                % outcomes as a function of z'
                N_next = (1-parms.s(iz))*parms.grid_N(iN) ...
                    + gesol.f(iN,ix,iz)*(1 - parms.grid_N(iN));
                N_next_notrunc = N_next;
                if N_next<Nlb
                    N_next = Nlb;
                elseif N_next>Nub
                    N_next = Nub;
                end
                x_next(:) = (1-parms.rho_x)*parms.x_bar + parms.rho_x*parms.grid_x(ix) ...
                    + parms.sigma_x(iN,iz)*(parms.grid_z(:) - expected_znext)/std_znext;
                
                phi_next(:) = phi_of_x(x_next);
                x_next(x_next<xlb) = xlb;
                x_next(x_next>xub) = xub;
                
                for iz_next = 1:Nz
                    C_next(iz_next) = interp2(N_mesh,x_mesh,gesol.C(:,:,iz_next)',N_next,x_next(iz_next));
                    S_next(iz_next) = interp2(N_mesh,x_mesh,gesol.S(:,:,iz_next)',N_next,x_next(iz_next));
                end
                
                % the SDF need not be truncated
                SDF(:) = parms.beta*((C_next./gesol.C(iN,ix,iz)).^(-parms.gamma))...
                    .*(zeta_fun(N_next_notrunc(:),phi_next(:))./zeta_fun(parms.grid_N(iN),phi_of_x(parms.grid_x(ix))));
                
%                 if any(~isfinite(SDF(:)))
%                     fprintf('got here.\n');
%                 end
                
                CV(iN,ix,iz) = parms.StateTransitionProbs(iz,:)*(SDF(:).*S_next(:));
                S(iN,ix,iz) = max(exp(parms.grid_z(iz)) - parms.b ...
                    + (1-parms.s(iz)-parms.eta*gesol.f(iN,ix,iz))*CV(iN,ix,iz),0);
                if ~isfinite(CV(iN,ix,iz))
                    fprintf('got here.\n');
                end
                if S(iN,ix,iz)==0
                    fprintf('got here.\n');
                end
                SDF_full(iN,ix,iz,:) = SDF;
                
            end
        end
    end
    
    err_CV = max(abs(gesol.CV(:) - CV(:)))
    err_S = max(abs(gesol.S(:) - S(:)))
    err_SDF = max(abs(gesol.spd1period(:) - SDF_full(:)))
    
    f_fun = @(x) (1+x.^(-parms.iota)).^(-1/parms.iota);
    g_fun = @(x) (1+x.^(parms.iota)).^(-1/parms.iota);
    finv_fun = @(f) (f.^(-parms.iota) - 1).^(-1/parms.iota);
    ginv_fun = @(g) (g.^(-parms.iota) - 1).^(1/parms.iota); % Solves for theta given g
    
    if isnan(parms.phi_match)
        fmax = ones(NN,Nx,Nz);
        gmin = zeros(NN,Nx,Nz);
    else
        fmax = ones(NN,Nx,Nz);
        for iN = 1:NN
            for iz = 1:Nz
                fmax(iN,:,iz) = max(min(1-parms.phi_match-parms.s(iz)+parms.s(iz)/(1-parms.grid_N(iN)),1),0);
            end
        end
        thetamax = finv_fun(fmax);
        gmin = g_fun(thetamax);
    end
    
    g = min(max(parms.kappa./((1-parms.eta)*CV),gmin),1);
    theta = ginv_fun(g);
    f = f_fun(theta);
    
    err_g = max(abs(g(:) - gesol.g(:)))
    err_theta = max(abs(theta(:) - gesol.theta(:)))
    err_f = max(abs(f(:) - gesol.f(:)))
    
    C = zeros(NN,Nx,iz);
    V = zeros(NN,Nx,iz);
    for iN = 1:NN
        for iz = 1:Nz
            V(iN,:,iz) = theta(iN,:,iz)*(1-parms.grid_N(iN));
            C(iN,:,iz) = exp(parms.grid_z(iz))*parms.grid_N(iN) ...
                - parms.kappa*V(iN,:,iz);
        end
    end
    
    err_V = max(abs(V(:) - gesol.V(:)))
    err_C = max(abs(C(:) - gesol.C(:)))

end